0%

【算法导论】矩阵乘法

发表于 更新于
本文字数： 12k 阅读时长 ≈ 17 分钟

离过年都不到十天了，还要等到这周五才能回家，想想也一年没回家了。从寒假开始到现在，已经有二十来天，这期间把2014年总结中的寒假计划也大多数完成了：The Element Of Style的阅读，三门数学课《随机过程》、《工程优化》、《数值分析》的算法实现。回家过年期间肯定不会写博客了，今天一看，这个月只写了三篇，于是乎今天必须再写一篇来完成这个月的基本工作量。言归正传，这篇文章写写选修课《算法设计》作业题中的矩阵乘法的三种方法。

矩阵乘法

  • 传统方法

    • 理论公式

    • 算法实现

      1
      2
      3
      4
      5
      6
      7
      8
      9
      10
      11
      12
      13
      14
      15
      16
      17
      void TraditionalMethod(float A[][N],float B[][N],float C[][N])//传统方法，三重循环
      {
          for(int i=0;i<N;i++)
              for(int j=0;j<N;j++)
                  C[i][j]=0;//之所以每次调用都清零，是因为前面是循环调用，如果只调用一次就不需要
          for(int i=0;i<N;i++)
          {
              for(int j=0;j<N;j++)
              {
                  for(int k=0;k<N;k++)
                  {
                      C[i][j]=C[i][j]+A[i][k]*B[k][j];
                  }
              }
          }

      }

  • 分块相乘法

    • 理论公式

    • 算法实现

      1
      2
      3
      4
      5
      6
      7
      8
      9
      10
      11
      12
      13
      14
      15
      16
      17
      18
      19
      20
      21
      22
      23
      24
      25
      26
      27
      28
      29
      30
      31
      32
      33
      34
      35
      36
      37
      38
      39
      40
      41
      42
      43
      44
      45
      46
      void BlockMatrix()//分块矩阵计算
      {     
            for(int i=0;i<N/2;i++)              
               for(int j=0;j<N/2;j++)
                  {
                      A11[i][j]=A[i][j];
                      A12[i][j]=A[i][j+N/2];
                      A21[i][j]=A[i+N/2][j];
                      A22[i][j]=A[i+N/2][j+N/2];
                      B11[i][j]=B[i][j];
                      B12[i][j]=B[i][j+N/2];
                      B21[i][j]=B[i+N/2][j];
                      B22[i][j]=B[i+N/2][j+N/2];

                      C11[i][j]=0;
                      C12[i][j]=0;
                      C21[i][j]=0;
                      C22[i][j]=0;
                  }       //将矩阵A和B式分为四块

               MATRIX_Multiply(N/2,A11,B11, AA);
               MATRIX_Multiply(N/2,A12,B21, BB);
               MATRIX_ADD(N/2,AA,BB,C11); //矩阵加法函数X+Y—>Z

               MATRIX_Multiply(N/2,A11,B12, AA);
               MATRIX_Multiply(N/2,A12,B22, BB);
               MATRIX_ADD(N/2,AA,BB,C12); //矩阵加法函数X+Y—>Z

               MATRIX_Multiply(N/2,A21,B11, AA);
               MATRIX_Multiply(N/2,A22,B21, BB);
               MATRIX_ADD(N/2,AA,BB,C21); //矩阵加法函数X+Y—>Z

               MATRIX_Multiply(N/2,A21,B12, AA);
               MATRIX_Multiply(N/2,A22,B22, BB);
               MATRIX_ADD(N/2,AA,BB,C22); //矩阵加法函数X+Y—>Z

          for(int i=0;i<N/2;i++)//将上面计算得到的结果放入结果矩阵C中
              for(int j=0;j<N/2;j++)
              {
                  C[i][j]=C11[i][j];
                  C[i][j+N/2]=C12[i][j];
                  C[i+N/2][j]=C21[i][j];
                  C[i+N/2][j+N/2]=C22[i][j];
              }                                            //计算结果送回C[N][N]

          }

  • Strassen法

    • 理论公式

    • 算法实现

      1
      2
      3
      4
      5
      6
      7
      8
      9
      10
      11
      12
      13
      14
      15
      16
      17
      18
      19
      20
      21
      22
      23
      24
      25
      26
      27
      28
      29
      30
      31
      32
      33
      34
      35
      36
      37
      38
      39
      40
      41
      42
      43
      44
      45
      46
      47
      48
      49
      50
      51
      52
      53
      54
      55
      56
      57
      58
      59
      60
      61
      62
      63
      64
      65
      66
      67
      68
      69
      70
      71
      void STRASSEN()  //STRASSEN函数
      {
          int i,j;//,x;

          for(i=0;i<N/2;i++)              
              for(j=0;j<N/2;j++)
              {
                  A11[i][j]=A[i][j];
                  A12[i][j]=A[i][j+N/2];
                  A21[i][j]=A[i+N/2][j];
                  A22[i][j]=A[i+N/2][j+N/2];
                  B11[i][j]=B[i][j];
                  B12[i][j]=B[i][j+N/2];
                  B21[i][j]=B[i+N/2][j];
                  B22[i][j]=B[i+N/2][j+N/2];
              }       //将矩阵A和B式分为四块




          MATRIX_SUB(N/2,B12,B22,BB);         
          MATRIX_Multiply(N/2,A11,BB,M1);

          MATRIX_ADD(N/2,A11,A12,AA);
          MATRIX_Multiply(N/2,AA,B22,M2);//M2=(A11+A12)B22

          MATRIX_ADD(N/2,A21,A22,AA);
          MATRIX_Multiply(N/2,AA,B11,M3);//M3=(A21+A22)B11

          MATRIX_SUB(N/2,B21,B11,BB);
          MATRIX_Multiply(N/2,A22,BB,M4);//M4=A22(B21-B11)

          MATRIX_ADD(N/2,A11,A22,AA);
          MATRIX_ADD(N/2,B11,B22,BB);
          MATRIX_Multiply(N/2,AA,BB,M5);//M5=(A11+A22)(B11+B22)


          MATRIX_SUB(N/2,A12,A22,AA);
          MATRIX_ADD(N/2,B21,B22,BB);
          MATRIX_Multiply(N/2,AA,BB,M6);//M6=(A12-A22)(B21+B22)

          MATRIX_SUB(N/2,A11,A21,AA);
          MATRIX_ADD(N/2,B11,B12,BB);
          MATRIX_Multiply(N/2,AA,BB,M7);//M7=(A11-A21)(B11+B12)
          //计算M1,M2,M3,M4,M5,M6,M7（递归部分）


          MATRIX_ADD(N/2,M5,M4,MM1);                
          MATRIX_SUB(N/2,M2,M6,MM2);
          MATRIX_SUB(N/2,MM1,MM2,C11);//C11=M5+M4-M2+M6

          MATRIX_ADD(N/2,M1,M2,C12);//C12=M1+M2

          MATRIX_ADD(N/2,M3,M4,C21);//C21=M3+M4

          MATRIX_ADD(N/2,M5,M1,MM1);
          MATRIX_ADD(N/2,M3,M7,MM2);
          MATRIX_SUB(N/2,MM1,MM2,C22);//C22=M5+M1-M3-M7

          for(i=0;i<N/2;i++)
              for(j=0;j<N/2;j++)
              {
                  C[i][j]=C11[i][j];
                  C[i][j+N/2]=C12[i][j];
                  C[i+N/2][j]=C21[i][j];
                  C[i+N/2][j+N/2]=C22[i][j];
              }                                            //计算结果送回C[N][N]



      }
  • 完整程序实现
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
#include <iostream>
#include<ctime>

using namespace std;

const int N=32; //常量N用来定义方阵的大小
void output(int n,float C[][N]); //函数声明部分
void TraditionalMethod(float A[][N],float B[][N],float C[][N]);//传统的矩阵相乘
void BlockMatrix();
void STRASSEN();
void MATRIX_Multiply(int n,float A[][N/2],float B[][N/2],float C[][N/2]);

float A[N][N];
float B[N][N];
float C[N][N];  //定义三个矩阵A,B,C

float A11[N/2][N/2],A12[N/2][N/2],A21[N/2][N/2],A22[N/2][N/2];
float B11[N/2][N/2],B12[N/2][N/2],B21[N/2][N/2],B22[N/2][N/2];
float C11[N/2][N/2],C12[N/2][N/2],C21[N/2][N/2],C22[N/2][N/2];
float M1[N/2][N/2],M2[N/2][N/2],M3[N/2][N/2],M4[N/2][N/2],M5[N/2][N/2],M6[N/2][N/2],M7[N/2][N/2];
float AA[N/2][N/2],BB[N/2][N/2],MM1[N/2][N/2],MM2[N/2][N/2];
void MATRIX_ADD(int n,float X[][N/2],float Y[][N/2],float Z[][N/2]); //矩阵加法函数X+Y—>Z

void main()
{
  //初始化，使相乘的两个矩阵都为全1矩阵
    for(int i=0;i<N;i++)
        for(int j=0;j<N;j++)
        {
            A[i][j]=1;
            B[i][j]=1;
        }
    //将结果矩阵C初始化为全0矩阵
    for(int i=0;i<N;i++)
        for(int j=0;j<N;j++)
            C[i][j]=0;

    clock_t start, finish; //用于计时
    double   duration; 
    int loop=0;

    cout<<"当矩阵为"<<N<<"×"<<N<<"，循环次数为10000时："<<endl<<endl;

    //--------使用传统方法--------------//
    start = clock();  
    while(loop<10000)//循环10000次,这里可以更改
    {
        loop++;
        TraditionalMethod(A,B,C);   //传统方法计算

    }
    finish = clock(); 
    duration = (double)(finish - start) / CLOCKS_PER_SEC; 
    cout<<"使用传统方法："<<endl;
    cout<<"所需时间为："<<duration<<endl<<endl;

  //  output(N,C);  //输出计算结果

    //---------使用分块矩阵乘法------------//
    start = clock();  
    loop=0;
    while(loop<10000)//循环10000次,这里可以更改
    {
        loop++;
        BlockMatrix();   //分块矩阵计算

    }
    finish = clock(); 
    duration = (double)(finish - start) / CLOCKS_PER_SEC; 
    cout<<"使用分块相乘方法："<<endl;
    cout<<"所需时间为："<<duration<<endl<<endl;
    //  output(N,C);  //输出计算结果

    //-------使用strassen方法-------------------//
    start = clock();  
    loop=0;
    while(loop<10000)//当时间非常小时，需要加大循环次数,这里可以更改
    {
        loop++;
        STRASSEN();   //调用STRASSEN函数计算      
    }
    finish = clock(); 
    duration = (double)(finish - start) / CLOCKS_PER_SEC; 
    cout<<"使用strassen方法："<<endl;
    cout<<"所需时间为："<<duration<<endl;
//  output(N,C);  //输出计算结果
}

void TraditionalMethod(float A[][N],float B[][N],float C[][N])//传统方法，三重循环
{
    for(int i=0;i<N;i++)
        for(int j=0;j<N;j++)
            C[i][j]=0;//之所以每次调用都清零，是因为前面是循环调用，如果只调用一次就不需要
    for(int i=0;i<N;i++)
    {
        for(int j=0;j<N;j++)
        {
            for(int k=0;k<N;k++)
            {
                C[i][j]=C[i][j]+A[i][k]*B[k][j];
            }
        }
    }

}


void BlockMatrix()//分块矩阵计算
{     
      for(int i=0;i<N/2;i++)              
         for(int j=0;j<N/2;j++)
            {
                A11[i][j]=A[i][j];
                A12[i][j]=A[i][j+N/2];
                A21[i][j]=A[i+N/2][j];
                A22[i][j]=A[i+N/2][j+N/2];
                B11[i][j]=B[i][j];
                B12[i][j]=B[i][j+N/2];
                B21[i][j]=B[i+N/2][j];
                B22[i][j]=B[i+N/2][j+N/2];

                C11[i][j]=0;
                C12[i][j]=0;
                C21[i][j]=0;
                C22[i][j]=0;
            }       //将矩阵A和B式分为四块

         MATRIX_Multiply(N/2,A11,B11, AA);
         MATRIX_Multiply(N/2,A12,B21, BB);
         MATRIX_ADD(N/2,AA,BB,C11); //矩阵加法函数X+Y—>Z

         MATRIX_Multiply(N/2,A11,B12, AA);
         MATRIX_Multiply(N/2,A12,B22, BB);
         MATRIX_ADD(N/2,AA,BB,C12); //矩阵加法函数X+Y—>Z

         MATRIX_Multiply(N/2,A21,B11, AA);
         MATRIX_Multiply(N/2,A22,B21, BB);
         MATRIX_ADD(N/2,AA,BB,C21); //矩阵加法函数X+Y—>Z

         MATRIX_Multiply(N/2,A21,B12, AA);
         MATRIX_Multiply(N/2,A22,B22, BB);
         MATRIX_ADD(N/2,AA,BB,C22); //矩阵加法函数X+Y—>Z

    for(int i=0;i<N/2;i++)//将上面计算得到的结果放入结果矩阵C中
        for(int j=0;j<N/2;j++)
        {
            C[i][j]=C11[i][j];
            C[i][j+N/2]=C12[i][j];
            C[i+N/2][j]=C21[i][j];
            C[i+N/2][j+N/2]=C22[i][j];
        }                                            //计算结果送回C[N][N]

    }



void output(int n,float C[][N]) //矩阵输出函数
{
    int i,j;
    cout<<"输出矩阵:"<<endl;
    for(i=0;i<n;i++)
    {
        cout<<endl;
        for(j=0;j<n;j++)
            cout<<C[i][j]<<"  ";
    }
    cout<<endl<<endl;

}


void MATRIX_Multiply(int n,float A[][N/2],float B[][N/2],float C[][N/2])//矩阵加法函数X*Y—>C
{


    for(int i=0;i<n;i++)
    {
        for(int j=0;j<n;j++)
        {
            C[i][j]=0;
            for(int k=0;k<n;k++)
            {
                C[i][j]=C[i][j]+A[i][k]*B[k][j];
            }
        }
    }

}



void MATRIX_ADD(int n,float X[][N/2],float Y[][N/2],float Z[][N/2]) //矩阵加法函数X+Y—>Z
{
    int i,j;
    for(i=0;i<n;i++)
        for(j=0;j<n;j++)
            Z[i][j]=X[i][j]+Y[i][j];
}

void MATRIX_SUB(int n,float X[][N/2],float Y[][N/2],float Z[][N/2]) //矩阵减法函数X-Y—>Z
{

    int i,j;
    for(i=0;i<n;i++)
        for(j=0;j<n;j++)
            Z[i][j]=X[i][j]-Y[i][j];

}

void STRASSEN()  //STRASSEN函数
{


    int i,j;//,x;

    for(i=0;i<N/2;i++)              
        for(j=0;j<N/2;j++)
        {
            A11[i][j]=A[i][j];
            A12[i][j]=A[i][j+N/2];
            A21[i][j]=A[i+N/2][j];
            A22[i][j]=A[i+N/2][j+N/2];
            B11[i][j]=B[i][j];
            B12[i][j]=B[i][j+N/2];
            B21[i][j]=B[i+N/2][j];
            B22[i][j]=B[i+N/2][j+N/2];
        }       //将矩阵A和B式分为四块




    MATRIX_SUB(N/2,B12,B22,BB);         
    MATRIX_Multiply(N/2,A11,BB,M1);

    MATRIX_ADD(N/2,A11,A12,AA);
    MATRIX_Multiply(N/2,AA,B22,M2);//M2=(A11+A12)B22

    MATRIX_ADD(N/2,A21,A22,AA);
    MATRIX_Multiply(N/2,AA,B11,M3);//M3=(A21+A22)B11

    MATRIX_SUB(N/2,B21,B11,BB);
    MATRIX_Multiply(N/2,A22,BB,M4);//M4=A22(B21-B11)

    MATRIX_ADD(N/2,A11,A22,AA);
    MATRIX_ADD(N/2,B11,B22,BB);
    MATRIX_Multiply(N/2,AA,BB,M5);//M5=(A11+A22)(B11+B22)


    MATRIX_SUB(N/2,A12,A22,AA);
    MATRIX_ADD(N/2,B21,B22,BB);
    MATRIX_Multiply(N/2,AA,BB,M6);//M6=(A12-A22)(B21+B22)

    MATRIX_SUB(N/2,A11,A21,AA);
    MATRIX_ADD(N/2,B11,B12,BB);
    MATRIX_Multiply(N/2,AA,BB,M7);//M7=(A11-A21)(B11+B12)
    //计算M1,M2,M3,M4,M5,M6,M7（递归部分）


    MATRIX_ADD(N/2,M5,M4,MM1);                
    MATRIX_SUB(N/2,M2,M6,MM2);
    MATRIX_SUB(N/2,MM1,MM2,C11);//C11=M5+M4-M2+M6

    MATRIX_ADD(N/2,M1,M2,C12);//C12=M1+M2

    MATRIX_ADD(N/2,M3,M4,C21);//C21=M3+M4

    MATRIX_ADD(N/2,M5,M1,MM1);
    MATRIX_ADD(N/2,M3,M7,MM2);
    MATRIX_SUB(N/2,MM1,MM2,C22);//C22=M5+M1-M3-M7

    for(i=0;i<N/2;i++)
        for(j=0;j<N/2;j++)
        {
            C[i][j]=C11[i][j];
            C[i][j+N/2]=C12[i][j];
            C[i+N/2][j]=C21[i][j];
            C[i+N/2][j+N/2]=C22[i][j];
        }                                            //计算结果送回C[N][N]   
}

运行结果如下图：

​ 图1

坚持原创技术分享，您的支持将鼓励我继续创作！